360 lines
14 KiB
Plaintext
360 lines
14 KiB
Plaintext
---
|
||
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"))
|
||
|
||
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"))
|
||
|
||
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, col = field)) +
|
||
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")
|
||
```
|
||
|