--- params: ref: "word-styles-reference-var1.docx" output_file: CI_report.docx report_date: "2024-08-28" data_dir: "Chemba" mail_day: "Wednesday" borders: TRUE output: # html_document: # toc: yes # df_print: paged word_document: reference_docx: !expr file.path("word-styles-reference-var1.docx") toc: yes editor_options: chunk_output_type: console --- ```{r setup, include=FALSE} #set de filename van de output # knitr::opts_chunk$set(echo = TRUE) # output_file <- params$output_file report_date <- params$report_date mail_day <- params$mail_day borders <- params$borders # # # # Activeer de renv omgeving # renv::activate() # renv::deactivate() # Optioneel: Herstel de omgeving als dat nodig is # Je kunt dit commentaar geven als je het normaal niet wilt uitvoeren # renv::restore() ``` ```{r libraries, message=FALSE, warning=FALSE, include=FALSE} knitr::opts_chunk$set(warning = FALSE, message = FALSE) 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) source("report_utils.R") # source(here("r_app", "report_utils.R")) ``` ```{r directories, message=FALSE, warning=FALSE, include=FALSE} project_dir <- params$data_dir source(here("r_app", "parameters_project.R")) log_message("Starting the R Markdown script") log_message(paste("mail_day params:", params$mail_day)) log_message(paste("report_date params:", params$report_date)) log_message(paste("mail_day variable:", mail_day)) # s2_dir <- "C:/Users/timon/Resilience BV/4002 CMD App - General/4002 CMD Team/4002 TechnicalData/04 WP2 technical/python/chemba_S2/" ``` ```{r week, message=FALSE, warning=FALSE, include=FALSE} Sys.setlocale("LC_TIME", "C") today <- as.character(report_date) mail_day_as_character <- as.character(mail_day) report_date_as_week_day <- weekdays(ymd(today)) days_of_week <- c("Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday") #als de index of report_date_as_week_day groter dan de index van de mail_day dan moet de week + 1 week <- week(today) log_message(paste("week", week, "today", today)) today_minus_1 <- as.character(ymd(today) - 7) today_minus_2 <- as.character(ymd(today) - 14) today_minus_3 <- as.character(ymd(today) - 21) log_message(paste("report_date_as_week_day", report_date_as_week_day)) log_message(paste("which(days_of_week == report_date_as_week_day)", which(days_of_week == report_date_as_week_day))) log_message(paste("mail_day_as_character", mail_day_as_character)) log_message(paste(" which(days_of_week == mail_day_as_character)", which(days_of_week == mail_day_as_character))) if (which(days_of_week == report_date_as_week_day) > which(days_of_week == mail_day_as_character)){ log_message("adjusting weeks because of mail day") week <- week(today) + 1 today_minus_1 <- as.character(ymd(today)) today_minus_2 <- as.character(ymd(today) - 7) today_minus_3 <- as.character(ymd(today) - 14) } # week <- week(today) # week <- 25 # today = "2024-06-21" #today = as.character(Sys.Date()) #week = lubridate::week(Sys.time()) ## week = 26 #title_var <- paste0("CI dashboard week ", week, " - all pivots dashboard using 3x3 meter resolution") subtitle_var <- paste("Report generated on", Sys.Date()) week_minus_1 <- week -1 # sprintf("%02d", week(today_minus_1)) week_minus_2 <- week -2 # sprintf("%02d", week(today_minus_2)) week_minus_3 <- week -3 # sprintf("%02d", week(today_minus_3)) week <- sprintf("%02d", week) year = year(today) year_1 = year(today_minus_1) year_2 = year(today_minus_2) year_3 = year(today_minus_3) ``` `r subtitle_var` \pagebreak # Explanation of the maps This PDF-dashboard shows the status of your fields on a weekly basis. We will show this in different ways: 1) The first way is with a general overview of field heterogeneity using ‘variation’ – a higher number indicates more differences between plants in the same field. 2) The second map shows a normal image of the latest week in colour, of the demo fields. 3) Then come the maps per field, which show the status for three weeks ago, two weeks ago, last week, and this week, as well as a percentage difference map between last week and this week. The percentage difference maps shows the relative difference in growth over the last week, with positive numbers showing growth, and negative numbers showing decline. 4) Below the maps are graphs that show how each pivot quadrant is doing, measured through the chlorophyll index. ```{r data, message=TRUE, warning=TRUE, include=FALSE} # get latest CI index # remove_pivots <- c("1.1", "1.12", "1.8", "1.9", "1.11", "1.14") CI_quadrant <- readRDS(here(cumulative_CI_vals_dir,"All_pivots_Cumulative_CI_quadrant_year_v2.rds"))# %>% # rename(pivot_quadrant = field) # path_to_week_current = here(weekly_CI_mosaic, paste0("week_",week, "_", year, ".tif")) # path_to_week_minus_1 = here(weekly_CI_mosaic, paste0("week_",week_minus_1, "_", year_1, ".tif")) # path_to_week_minus_2 = here(weekly_CI_mosaic, paste0("week_",week_minus_2, "_", year_2, ".tif")) # path_to_week_minus_3 = here(weekly_CI_mosaic, paste0("week_",week_minus_3, "_", year_3, ".tif")) path_to_week_current = get_week_path(weekly_CI_mosaic, today, 0) path_to_week_minus_1 = get_week_path(weekly_CI_mosaic, today, -1) path_to_week_minus_2 = get_week_path(weekly_CI_mosaic, today, -2) path_to_week_minus_3 = get_week_path(weekly_CI_mosaic, today, -3) log_message("required mosaic paths") log_message(paste("path to week current",path_to_week_current)) log_message(paste("path to week minus 1",path_to_week_minus_1)) log_message(paste("path to week minus 2",path_to_week_minus_2)) log_message(paste("path to week minus 3",path_to_week_minus_3)) CI <- brick(path_to_week_current) %>% subset("CI") CI_m1 <- brick(path_to_week_minus_1) %>% subset("CI") CI_m2 <- brick(path_to_week_minus_2) %>% subset("CI") CI_m3 <- brick(path_to_week_minus_3) %>% subset("CI") # last_week_dif_raster <- ((CI - CI_m1) / CI_m1) * 100 last_week_dif_raster_abs <- (CI - CI_m1) ``` ```{r data_129, message=TRUE, warning=TRUE, include=FALSE} three_week_dif_raster_abs <- (CI - CI_m3) ``` ```{r data_132, message=TRUE, warning=TRUE, include=FALSE} # AllPivots0 <-st_read(here(data_dir_project, "pivot.geojson")) # AllPivots0$pivot <- factor(AllPivots0$pivot, levels = c("1.1", "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" ,"2.1", "2.2", "2.3" , "2.4", "2.5", "3.1", "3.2", "3.3", "4.1", "4.2", "4.3", "4.4", "4.5", "4.6", "5.1" ,"5.2", "5.3", "5.4", "6.1", "6.2", "DL1.1", "DL1.3")) AllPivots0 <- field_boundaries_sf # joined_spans <-st_read(here(data_dir_project, "span.geojson")) %>% st_transform(crs(AllPivots0)) # pivots_dates <- readRDS(here(harvest_dir, "harvest_data_new")) # pivots_dates$pivot <- factor(pivots_dates$pivot, levels = c("1.1", "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" ,"2.1", "2.2", "2.3" , "2.4", "2.5", "3.1", "3.2", "3.3", "4.1", "4.2", "4.3", "4.4", "4.5", "4.6", "5.1" ,"5.2", "5.3", "5.4", "6.1", "6.2", "DL1.1", "DL1.3")) #AllPivots <- merge(AllPivots0, harvesting_data, by = c("field", "sub_field")) #%>% #rename(field = pivot, sub_field = pivot_quadrant) #%>% select(-pivot.y) #head(AllPivots) #AllPivots_merged <- AllPivots %>% #dplyr::select(field, sub_field, sub_area) %>% unique() %>% # group_by(field) %>% summarise(sub_area = first(sub_area)) #AllPivots_merged <- st_transform(AllPivots_merged, crs = proj4string(CI)) #pivot_names <- unique(CI_quadrant$field) ``` ```{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 RGB_raster <- brick(RGB_raster) # RGB_raster <- brick(here("planet", paste0("week_",week, ".tif"))) tm_shape(RGB_raster, unit = "m") + tm_rgb(r=1, g=2, b=3, max.value = 255) + tm_layout(main.title = paste0("RGB image of the fields - week ", week), main.title.position = 'center') + tm_scale_bar(position = c("right", "top"), text.color = "white") + tm_compass(position = c("right", "top"), text.color = "white") + tm_shape(AllPivots0)+ tm_borders( col = "white") + tm_text("pivot_quadrant", size = .6, col = "white") ``` \newpage ```{r ci_overzicht_kaart, echo=FALSE, fig.height=7.3, fig.width=9, message=FALSE, warning=FALSE} tm_shape(CI, unit = "m")+ tm_raster(breaks = c(0,0.5,1,2,3,4,5,6,7,Inf), palette = "RdYlGn", midpoint = NA,legend.is.portrait = F) + tm_layout(legend.outside = TRUE,legend.outside.position = "bottom",legend.show = T, main.title = "Overview all fields (CI)")+ tm_scale_bar(position = c("right", "top"), text.color = "black") + tm_compass(position = c("right", "top"), text.color = "black") + tm_shape(AllPivots0)+ tm_borders( col = "black") + tm_text("sub_field", size = .6, col = "black") ``` \newpage ```{r ci_diff_kaart, echo=FALSE, fig.height=7.3, fig.width=9, message=FALSE, warning=FALSE} tm_shape(last_week_dif_raster_abs, unit = "m")+ tm_raster(breaks = c(-3,-2,-1,0,1,2, 3), palette = "RdYlGn", midpoint = NA,legend.is.portrait = F) + tm_layout(legend.outside = TRUE,legend.outside.position = "bottom",legend.show = T, main.title = "Overview all fields - CI difference")+ tm_scale_bar(position = c("right", "top"), text.color = "black") + tm_compass(position = c("right", "top"), text.color = "black") + tm_shape(AllPivots0)+ tm_borders( col = "black") + tm_text("sub_field", size = .6, col = "black") ``` \newpage ```{r plots_ci_estate, eval=TRUE, fig.height=3.8, fig.width=10, message=FALSE,echo=FALSE, warning=FALSE, include=TRUE, results='asis'} AllPivots_merged <- AllPivots0 %>% dplyr::group_by(field) %>% summarise() walk(AllPivots_merged$field, ~ { cat("\n") # Add an empty line for better spacing ci_plot(.x) cat("\n") cum_ci_plot(.x) }) ``` ```{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 # Iterate over each subgroup for (subgroup in unique(pivots_grouped$sub_area)) { cat("# HELLO!!!") print(" PRINT") # cat("\n") # cat("# Subgroup: ", subgroup, "\n") # Add a title for the subgroup subset_data <- filter(pivots_grouped, sub_area == subgroup) # cat("\\pagebreak") walk(subset_data$field, ~ { # cat("\n") # Add an empty line for better spacing ci_plot(.x) # cat("\n") cum_ci_plot(.x) # cat("\n") }) } ``` # Yield prediction The below table shows estimates of the biomass if you would harvest them now. ```{r message=FALSE, warning=FALSE, include=FALSE} #CI_quadrant <- readRDS(here(cumulative_CI_vals_dir,"All_pivots_Cumulative_CI_quadrant_year_v2.rds")) CI_quadrant <- readRDS(here(cumulative_CI_vals_dir,"All_pivots_Cumulative_CI_quadrant_year_v2.rds")) %>% group_by(model) %>% # Group by model or other relevant columns fill(field, sub_field, .direction = "downup") %>% # Fill down then up within each group ungroup() harvesting_data <- harvesting_data %>% rename(season = year) CI_and_yield <- left_join(CI_quadrant , harvesting_data, by = c("field", "sub_field", "season")) %>% #filter(!is.na(tonnage_ha)) %>% group_by(sub_field, season) %>% slice(which.max(DOY)) %>% dplyr::select(field, sub_field, tonnage_ha, cumulative_CI, DOY, season, sub_area) %>% mutate(CI_per_day = cumulative_CI/ DOY) predictors <- c( "cumulative_CI" , "DOY" ,"CI_per_day" ) response <- "tonnage_ha" # CI_and_yield_test <- as.data.frame(CI_and_yield_test) CI_and_yield_test <- CI_and_yield %>% as.data.frame() %>% filter(!is.na(tonnage_ha)) CI_and_yield_validation <- CI_and_yield_test prediction_yields <- CI_and_yield %>% as.data.frame() %>% filter(is.na(tonnage_ha)) ctrl <- trainControl(method="cv", savePredictions = TRUE, allowParallel= TRUE, number = 5, verboseIter = TRUE) set.seed(202) model_ffs_rf <- ffs( CI_and_yield_test[,predictors], CI_and_yield_test[,response], method="rf" , trControl=ctrl, importance=TRUE, withinSE = TRUE, tuneLength = 5, na.rm = TRUE ) # Function to prepare predictions prepare_predictions <- function(predictions, newdata) { return(predictions %>% as.data.frame() %>% rename(predicted_Tcha = ".") %>% mutate(sub_field = newdata$sub_field, field = newdata$field, Age_days = newdata$DOY, total_CI = round(newdata$cumulative_CI, 0), predicted_Tcha = round(predicted_Tcha, 0), season = newdata$season) %>% dplyr::select(field, sub_field, Age_days, total_CI, predicted_Tcha, season) %>% left_join(., newdata, by = c("field", "sub_field", "season"))) } # Predict yields for the validation dataset pred_ffs_rf <- prepare_predictions(predict(model_ffs_rf, newdata = CI_and_yield_validation), CI_and_yield_validation) # Predict yields for the current season pred_rf_current_season <- prepare_predictions(predict(model_ffs_rf, newdata = prediction_yields), prediction_yields) %>% filter(Age_days > 300) %>% mutate(CI_per_day = round(total_CI / Age_days, 1)) ``` ```{r yield_plaatjes, echo=FALSE } ggplot(pred_ffs_rf, aes(y = predicted_Tcha, x = tonnage_ha)) + geom_point(size = 2, alpha = 0.6) + # Adjust point size and transparency geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "red") + # Reference line scale_x_continuous(limits = c(0, 200)) + scale_y_continuous(limits = c(0, 200)) + labs(title = "Model Performance: \nPredicted vs Actual Tonnage/ha", x = "Actual tonnage/ha (Tcha)", y = "Predicted tonnage/ha (Tcha)") + theme_minimal() ggplot(pred_rf_current_season, aes(x = Age_days, y = predicted_Tcha)) + geom_point(size = 2, alpha = 0.6) + # Adjust point size and transparency labs(title = "Predicted Yields for Fields Over 300 Days \nOld Yet to Be Harvested", x = "Age (days)", y = "Predicted tonnage/ha (Tcha)") + # facet_wrap(~sub_area) + scale_y_continuous(limits = c(0, 200)) + # Optional: Set limits for y-axis theme_minimal() knitr::kable(pred_rf_current_season, digits = 0, caption = "Predicted Tonnage/ha for Fields Over 300 Days Old") ```