From d0ff25acf926add66232927094bdcf44ff3e3849 Mon Sep 17 00:00:00 2001 From: Timon Date: Mon, 16 Dec 2024 16:26:34 +0100 Subject: [PATCH] timely save --- r_app/plot_testing.R | 117 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 117 insertions(+) create mode 100644 r_app/plot_testing.R diff --git a/r_app/plot_testing.R b/r_app/plot_testing.R new file mode 100644 index 0000000..3dc5b97 --- /dev/null +++ b/r_app/plot_testing.R @@ -0,0 +1,117 @@ +library(here) +library(sf) +library(tidyverse) +library(tmap) +library(lubridate) +library(exactextractr) +library(zoo) +library(raster) +library(terra) +library(rsample) +library(caret) +library(randomForest) +library(CAST) + +project_dir <- "chemba" +source(here("r_app", "parameters_project.R")) +cumulative_CI_vals_dir <- "C:/Users/timon/Resilience BV/4020 SCane ESA DEMO - Documenten/General/4020 SCDEMO Team/4020 TechnicalData/WP3/smartcane/laravel_app/storage/app/chemba/Data/extracted_ci/cumulative_vals" +CI_quadrant <- readRDS(here(cumulative_CI_vals_dir,"All_pivots_Cumulative_CI_quadrant_year_v2.rds"))# %>% + +cum_ci_plot <- function(pivotName, plot_type = "value", facet_on = TRUE, x_axis = "days") { + + # pivotName = "3a13" + data_ci <- CI_quadrant %>% filter(field == pivotName) + + if (nrow(data_ci) == 0) { + return(cum_ci_plot2(pivotName)) # Return an empty data frame if no data is found + } + + data_ci2 <- data_ci %>% + mutate(CI_rate = cumulative_CI / DOY, + week = week(Date), + week_from_doy = DOY / 7) %>% + group_by(field) %>% + mutate(mean_CIrate_rolling_10_days = rollapplyr(CI_rate, width = 10, FUN = mean, partial = TRUE), + mean_rolling_10_days = rollapplyr(value, width = 10, FUN = mean, partial = TRUE)) + + data_ci2 <- data_ci2 %>% mutate(season = as.factor(season)) + + date_preperation_perfect_pivot <- data_ci2 %>% + group_by(season) %>% + summarise(min_date = min(Date), + max_date = max(Date), + days = max_date - min_date) + + unique_seasons <- sort(unique(date_preperation_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") + + if (facet_on) { + g <- ggplot(data = data_ci2 %>% filter(season %in% unique_seasons)) + + facet_wrap(~season, scales = "free_x") + + geom_line(aes_string(x = "Date", y = y_aesthetic, col = "sub_field", group = "sub_field")) + + labs(title = paste("Plot of", y_label, "for Pivot", pivotName), + color = "Field Name", + y = y_label) + + scale_x_date(date_breaks = "1 month", date_labels = "%m-%Y") + + 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)) + + guides(color = guide_legend(nrow = 2, byrow = TRUE)) + } else if (x_axis == "weeks") { + g <- ggplot(data = data_ci2 %>% filter(season %in% unique_seasons)) + + facet_wrap(~sub_field, nrow=1) + + geom_line(aes_string(x = "week_from_doy", y = y_aesthetic, col = "season", group = "season")) + + labs(title = paste("Plot of", y_label, "for Pivot", pivotName), + color = "Season", + y = y_label, + x = "Week of Year") + + 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)) + + guides(color = guide_legend(nrow = 2, byrow = TRUE)) + } else { + g <- ggplot(data = data_ci2 %>% filter(season %in% unique_seasons)) + + facet_wrap(~sub_field, nrow=1) + + geom_line(aes_string(x = "DOY", y = y_aesthetic, col = "season", group = "season")) + + labs(title = paste("Plot of", y_label, "for Pivot", pivotName), + color = "Season", + y = y_label, + x = "Age of Crop (Days)") + + 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)) + + guides(color = guide_legend(nrow = 2, byrow = TRUE)) + } + + return(g) +} + +# facet_on FALSE +g1 <- cum_ci_plot("1.16", "value", FALSE, "days") +g2 <- cum_ci_plot("3a13", "cumulative_CI", FALSE, "days") +g3 <- cum_ci_plot("3a13", "CI_rate", FALSE, "days") +g7 <- cum_ci_plot("3a13", "value", FALSE, "weeks") +g8 <- cum_ci_plot("3a13", "cumulative_CI", FALSE, "weeks") +g9 <- cum_ci_plot("3a13", "CI_rate", FALSE, "weeks") + +# facet_on TRUE +g4 <- cum_ci_plot("3a13", "value", TRUE, "days") +g5 <- cum_ci_plot("3a13", "cumulative_CI", TRUE, "days") +g6 <- cum_ci_plot("3a13", "CI_rate", TRUE, "days") +