diff --git a/r_app/CI_report_dashboard_planet.Rmd b/r_app/CI_report_dashboard_planet.Rmd index 3ed7508..2982dfc 100644 --- a/r_app/CI_report_dashboard_planet.Rmd +++ b/r_app/CI_report_dashboard_planet.Rmd @@ -2,8 +2,8 @@ params: ref: "word-styles-reference-var1.docx" output_file: CI_report.docx - report_date: "2024-04-18" - data_dir: "Sony" + report_date: "2024-08-28" + data_dir: "Chemba" mail_day: "Wednesday" borders: TRUE output: @@ -50,6 +50,9 @@ library(caret) library(randomForest) library(CAST) +source(here("r_app/report_utils.R")) +# source("report_utils.R") + # Define the log file path log_file <- here("laravel_app/storage/app/rmd_log.txt") @@ -202,197 +205,6 @@ AllPivots0 <- field_boundaries_sf ``` -```{r eval=FALSE, include=FALSE} -pivot_stats_q <- cbind(AllPivots0, round(exact_extract(CI, AllPivots0, c("coefficient_of_variation", "mean"), default_value = -9999),2)) %>% - st_drop_geometry() %>% as_tibble() - -hetero_pivots0 <- merge(AllPivots, pivot_stats_q %>% dplyr::select(-hectares, -radius, -pivot), by = "pivot_quadrant") - -hetero_pivots <- hetero_pivots0 %>% #dplyr::filter(variable %in% "CV") %>% - mutate(class = case_when(coefficient_of_variation <= 0.001 ~ "Missing data", - coefficient_of_variation >= 0.002 & coefficient_of_variation < 0.1 ~ "Homogeneous", - coefficient_of_variation >= 0.1 & coefficient_of_variation < 0.2~ "Somewhat Homogeneous", - coefficient_of_variation >= 0.2 & coefficient_of_variation < 0.4 ~ "Somewhat Heterogeneous", - coefficient_of_variation >= 0.4 ~ "Heterogeneous")) %>% - mutate(class = as.factor(class)) -# hetero_pivots %>% filter(pivot == "1.3") -hetero_pivots$class <- factor(hetero_pivots$class, levels = c("Missing data", "Homogeneous","Somewhat Homogeneous", "Somewhat Heterogeneous", "Heterogeneous" )) -# hetero_pivots %>% select(pivot_quadrant, class, Age, coefficient_of_variation , mean) %>% view() - -# hetero_pivots %>% filter(class == "Somewhat Heterogeneous" ) %>% select(pivot_quadrant, class, Age, coefficient_of_variation , mean) - - - -Mypal <- c('#dcdcdc', '#008000','#8db600','#FFC300','#F22222') - -hetero_plot <- function(data){ - # map <- - tm_shape(data) + tm_polygons(col = "class", palette=Mypal) + - tm_text("pivot_quadrant", size = 1/2) + - tm_layout(main.title=paste0("Homogeneity of pivot quadrants, week ", week, " 2022"),main.title.position = "center")+ - tm_compass(position = c("center", "top")) + - tm_scale_bar(position = c("center", "top")) - - # print(map, width = 20, units = "cm") - -} - -``` -\newpage -```{r eval=FALSE, fig.height=7, fig.width=10, message=FALSE, warning=FALSE, include=FALSE} -hetero_plot(hetero_pivots) -``` - - -```{r function, message=FALSE, warning=FALSE, include=FALSE} - -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)) -} - -create_CI_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week, age, borders = FALSE){ - map <- tm_shape(pivot_raster, unit = "m") + - tm_raster(breaks = c(0,0.5,1,2,3,4,5,6,7,Inf), palette = "RdYlGn",legend.is.portrait = legend_is_portrait ,midpoint = NA) + - tm_layout(main.title = paste0("\nMax CI week ", week,"\n", age, " weeks old"), - main.title.size = 0.7, legend.show = show_legend) - - 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) -} - -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){ - map <- tm_shape(pivot_raster, unit = "m") + - tm_raster(breaks = c(-3,-2,-1,0,1,2,3), palette = "RdYlGn",legend.is.portrait = legend_is_portrait, midpoint = 0, title = "CI difference") + - tm_layout(main.title = paste0("CI change week ", week_1, " - week ", week_2, "\n", age, " weeks old"), - main.title.size = 0.7, legend.show = show_legend) - - 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) -} - -ci_plot <- function(pivotName){ - # pivotName = "1.1" - pivotShape <- AllPivots0 %>% terra::subset(field %in% pivotName) %>% st_transform(crs(CI)) - age <- harvesting_data %>% dplyr::filter(field %in% pivotName) %>% sort("year") %>% tail(., 1) %>% dplyr::select(age) %>% unique() %>% pull() %>% round() - - AllPivots2 <- AllPivots0 %>% dplyr::filter(field %in% pivotName) - - singlePivot <- CI %>% crop(., pivotShape) %>% mask(., pivotShape) - - singlePivot_m1 <- CI_m1 %>% crop(., pivotShape) %>% mask(., pivotShape) - singlePivot_m2 <- CI_m2 %>% crop(., pivotShape) %>% mask(., pivotShape) - # singlePivot_m3 <- CI_m3 %>% crop(., pivotShape) %>% mask(., pivotShape) - - abs_CI_last_week <- last_week_dif_raster_abs %>% crop(., pivotShape) %>% mask(., pivotShape) - abs_CI_three_week <- three_week_dif_raster_abs %>% crop(., pivotShape) %>% mask(., pivotShape) - - planting_date <- harvesting_data %>% dplyr::filter(field %in% pivotName) %>% ungroup() %>% dplyr::select(season_start) %>% unique() - - joined_spans2 <- AllPivots0 %>% st_transform(crs(pivotShape)) %>% dplyr::filter(field %in% pivotName) #%>% unique() %>% st_crop(., pivotShape) - - CImap_m2 <- create_CI_map(singlePivot_m2, AllPivots2, joined_spans2, show_legend= T, legend_is_portrait = T, week = week_minus_2, age = age -2, borders = borders) - CImap_m1 <- create_CI_map(singlePivot_m1, AllPivots2, joined_spans2, show_legend= F, legend_is_portrait = F, week = week_minus_1, age = age -1, borders = borders) - CImap <- create_CI_map(singlePivot, AllPivots2, joined_spans2, show_legend= F, legend_is_portrait = F, week = week, age = age, borders = borders) - - - CI_max_abs_last_week <- create_CI_diff_map(abs_CI_last_week,AllPivots2, joined_spans2, show_legend = T, legend_is_portrait = T, 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 = T, legend_is_portrait = T, week_1 = week, week_2 = week_minus_3, age = age, borders = borders) - - tst <- tmap_arrange(CImap_m2, CImap_m1, CImap,CI_max_abs_last_week, CI_max_abs_three_week, nrow = 1) - - cat(paste("## Field", pivotName, "-", age, "weeks after planting/harvest", "\n")) - # cat("\n") - # cat('

Pivot', pivotName, '- week', week, '-', age$Age, 'weeks after planting/harvest

') - # cat(paste("# Pivot",pivots$pivot[i],"\n")) - print(tst) - - } - -cum_ci_plot <- function(pivotName){ - - # pivotName = "1.1" - 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))%>% group_by(field) %>% - mutate(mean_rolling10 = rollapplyr(CI_rate , width = 10, FUN = mean, partial = TRUE)) - - 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 <- unique(date_preperation_perfect_pivot$season) - - g <- ggplot(data= data_ci2 %>% filter(season %in% unique_seasons)) + - facet_wrap(~season, scales = "free_x") + - geom_line( aes(Date, mean_rolling10, col = sub_field, group = sub_field)) + - labs(title = paste("14 day rolling MEAN CI rate - Pivot ", pivotName))+ - 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)) - - subchunkify(g, 3.2, 10) - -} - -cum_ci_plot2 <- function(pivotName){ - 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(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) -} - - -``` - ```{r eval=FALSE, fig.height=7.2, fig.width=10, message=FALSE, warning=FALSE, include=FALSE} RGB_raster <- list.files(paste0(s2_dir,week),full.names = T, pattern = ".tiff", recursive = TRUE)[1] #use pattern = '.tif$' or something else if you have multiple files in this folder @@ -439,13 +251,6 @@ tm_shape(CI, unit = "m")+ \newpage ```{r plots_ci_estate, eval=TRUE, fig.height=3.8, fig.width=10, message=FALSE,echo=FALSE, warning=FALSE, include=TRUE, results='asis'} -# # pivots <- AllPivots_merged %>% filter(pivot != c("1.1", "1.17")) -# pivots_estate <- AllPivots_merged %>% filter(field %in% c("6.2")) #, "1.2", "1.3", "1.4")) #, "1.6", "1.7", "1.8", "1.9", "1.10", "1.11", "1.12", "1.13", "1.14" , "1.16" , "1.17" , "1.18" , "6.1", "6.2", "DL1.1", "DL1.3")) %>% filter(pivot != "1.17") - - -# pivots <- AllPivots_merged %>% filter(pivot != c("1.1", "1.17")) -# pivots_estate <- AllPivots_merged %>% filter(pivot %in% c("1.1", "1.2", "1.7")) %>% filter(pivot != "1.17") - AllPivots_merged <- AllPivots0 %>% dplyr::group_by(field) %>% summarise() walk(AllPivots_merged$field, ~ { @@ -457,9 +262,7 @@ walk(AllPivots_merged$field, ~ { ``` ```{r looping_over_sub_area, echo=FALSE, fig.height=3.8, fig.width=10, message=FALSE, warning=FALSE, results='asis', eval=FALSE} -pivots_grouped <- AllPivots0 # %>% -# group_by(sub_area) %>% -# arrange(sub_area) # Optional: arrange the groups alphabetically by sub_area +pivots_grouped <- AllPivots0 # Iterate over each subgroup for (subgroup in unique(pivots_grouped$sub_area)) { diff --git a/r_app/interpolate_growth_model.R b/r_app/interpolate_growth_model.R index b5f43eb..1cf843d 100644 --- a/r_app/interpolate_growth_model.R +++ b/r_app/interpolate_growth_model.R @@ -17,7 +17,7 @@ project_dir <- as.character(args[1]) # Controleer of data_dir een geldige waarde is if (!is.character(project_dir)) { - project_dir <- "sony" + project_dir <- "chemba" } @@ -133,8 +133,6 @@ CI_all <- CI_all %>% group_by(model) %>% mutate(CI_per_day = FitData - lag(FitDa cumulative_CI = cumsum(FitData)) -CI_all <- CI_all %>% group_by(model) %>% mutate(CI_per_day = FitData - lag(FitData), - cumulative_CI = cumsum(FitData)) message('CI_all cumulative') head(CI_all) diff --git a/r_app/report_utils.R b/r_app/report_utils.R new file mode 100644 index 0000000..5aff7f2 --- /dev/null +++ b/r_app/report_utils.R @@ -0,0 +1,145 @@ +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)) +} + +create_CI_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week, age, borders = FALSE){ + map <- tm_shape(pivot_raster, unit = "m") + + tm_raster(breaks = c(0,0.5,1,2,3,4,5,6,7,Inf), palette = "RdYlGn",legend.is.portrait = legend_is_portrait ,midpoint = NA) + + tm_layout(main.title = paste0("\nMax CI week ", week,"\n", age, " weeks old"), + main.title.size = 0.7, legend.show = show_legend) + + 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) +} + +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){ + map <- tm_shape(pivot_raster, unit = "m") + + tm_raster(breaks = c(-3,-2,-1,0,1,2,3), palette = "RdYlGn",legend.is.portrait = legend_is_portrait, midpoint = 0, title = "CI difference") + + tm_layout(main.title = paste0("CI change week ", week_1, " - week ", week_2, "\n", age, " weeks old"), + main.title.size = 0.7, legend.show = show_legend) + + 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) +} + +ci_plot <- function(pivotName){ + # pivotName = "1.1" + pivotShape <- AllPivots0 %>% terra::subset(field %in% pivotName) %>% st_transform(crs(CI)) + age <- harvesting_data %>% dplyr::filter(field %in% pivotName) %>% sort("year") %>% tail(., 1) %>% dplyr::select(age) %>% unique() %>% pull() %>% round() + + AllPivots2 <- AllPivots0 %>% dplyr::filter(field %in% pivotName) + + singlePivot <- CI %>% crop(., pivotShape) %>% mask(., pivotShape) + + singlePivot_m1 <- CI_m1 %>% crop(., pivotShape) %>% mask(., pivotShape) + singlePivot_m2 <- CI_m2 %>% crop(., pivotShape) %>% mask(., pivotShape) + # singlePivot_m3 <- CI_m3 %>% crop(., pivotShape) %>% mask(., pivotShape) + + abs_CI_last_week <- last_week_dif_raster_abs %>% crop(., pivotShape) %>% mask(., pivotShape) + abs_CI_three_week <- three_week_dif_raster_abs %>% crop(., pivotShape) %>% mask(., pivotShape) + + planting_date <- harvesting_data %>% dplyr::filter(field %in% pivotName) %>% ungroup() %>% dplyr::select(season_start) %>% unique() + + joined_spans2 <- AllPivots0 %>% st_transform(crs(pivotShape)) %>% dplyr::filter(field %in% pivotName) #%>% unique() %>% st_crop(., pivotShape) + + CImap_m2 <- create_CI_map(singlePivot_m2, AllPivots2, joined_spans2, show_legend= T, legend_is_portrait = T, week = week_minus_2, age = age -2, borders = borders) + CImap_m1 <- create_CI_map(singlePivot_m1, AllPivots2, joined_spans2, show_legend= F, legend_is_portrait = F, week = week_minus_1, age = age -1, borders = borders) + CImap <- create_CI_map(singlePivot, AllPivots2, joined_spans2, show_legend= F, legend_is_portrait = F, week = week, age = age, borders = borders) + + + CI_max_abs_last_week <- create_CI_diff_map(abs_CI_last_week,AllPivots2, joined_spans2, show_legend = T, legend_is_portrait = T, 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 = T, legend_is_portrait = T, week_1 = week, week_2 = week_minus_3, age = age, borders = borders) + + tst <- tmap_arrange(CImap_m2, CImap_m1, CImap,CI_max_abs_last_week, CI_max_abs_three_week, nrow = 1) + + cat(paste("## Field", pivotName, "-", age, "weeks after planting/harvest", "\n")) + # cat("\n") + # cat('

Pivot', pivotName, '- week', week, '-', age$Age, 'weeks after planting/harvest

') + # cat(paste("# Pivot",pivots$pivot[i],"\n")) + print(tst) + +} + +cum_ci_plot <- function(pivotName){ + + # pivotName = "1.1" + 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))%>% group_by(field) %>% + mutate(mean_rolling10 = rollapplyr(CI_rate , width = 10, FUN = mean, partial = TRUE)) + + 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 <- unique(date_preperation_perfect_pivot$season) + + g <- ggplot(data= data_ci2 %>% filter(season %in% unique_seasons)) + + facet_wrap(~season, scales = "free_x") + + geom_line( aes(Date, mean_rolling10, col = sub_field, group = sub_field)) + + labs(title = paste("14 day rolling MEAN CI rate - Pivot ", pivotName), + color = "Field name")+ + 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)) + + subchunkify(g, 3.2, 10) + +} + +cum_ci_plot2 <- function(pivotName){ + 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(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) +} \ No newline at end of file