541 lines
24 KiB
Plaintext
541 lines
24 KiB
Plaintext
---
|
||
# title: paste0("CI report week ", week, " - all pivots from ", last_tuesday, " to ", today)
|
||
params:
|
||
ref: word-styles-reference-var1.docx
|
||
output_file: "CI_report.docx"
|
||
report_date: "2023-12-12"
|
||
output:
|
||
word_document:
|
||
reference_docx: !expr file.path("word-styles-reference-var1.docx")
|
||
# toc: true
|
||
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
|
||
|
||
|
||
# Activeer de renv omgeving
|
||
renv::activate()
|
||
|
||
# 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(rsample)
|
||
library(caret)
|
||
library(randomForest)
|
||
library(CAST)
|
||
```
|
||
|
||
```{r directories, message=FALSE, warning=FALSE, include=FALSE}
|
||
laravel_storage_dir <- here("../laravel_app/storage/app/chemba")
|
||
data_dir <- here(laravel_storage_dir, "../Data")
|
||
extracted_CI_dir <- here(data_dir, "extracted_ci")
|
||
daily_CI_vals_dir <- here(extracted_CI_dir, "daily_vals")
|
||
cumulative_CI_vals_dir <- here(extracted_CI_dir, "cumulative_vals")
|
||
harvest_dir <- here(data_dir, "HarvestData")
|
||
|
||
weekly_CI_mosaic <- here(laravel_storage_dir, "weekly_mosaic")
|
||
|
||
|
||
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}
|
||
# week <- 5
|
||
#today = "2023-12-12"
|
||
today <- as.character(report_date)
|
||
week <- week(today)
|
||
|
||
|
||
#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", 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)
|
||
|
||
week_minus_1 <- sprintf("%02d", week(today_minus_1))
|
||
week_minus_2 <- sprintf("%02d", week(today_minus_2))
|
||
week_minus_3 <- sprintf("%02d", week(today_minus_3))
|
||
|
||
year = year(today)
|
||
year_2 = year(today_minus_1)
|
||
year_3 = year(today_minus_2)
|
||
year_4 = 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)
|
||
|
||
message("STOP - check ci name in layer")
|
||
CI <- brick(here(weekly_CI_mosaic, paste0("week_",week, "_", year, ".tif"))) %>% subset("CI")
|
||
CI_m1 <- brick(here(weekly_CI_mosaic, paste0("week_",week_minus_1, "_", year_2, ".tif"))) %>% subset("CI")
|
||
CI_m2 <- brick(here(weekly_CI_mosaic, paste0("week_",week_minus_2, "_", year_3, ".tif"))) %>% subset("CI")
|
||
CI_m3 <- brick(here(weekly_CI_mosaic, paste0("week_",week_minus_3, "_", year_4, ".tif"))) %>% subset("CI")
|
||
|
||
# last_week_dif_raster <- ((CI - CI_m1) / CI_m1) * 100
|
||
last_week_dif_raster_abs <- (CI - CI_m1)
|
||
three_week_dif_raster_abs <- (CI - CI_m3)
|
||
|
||
AllPivots0 <-st_read(here(data_dir, "pivot_20210625.geojson"))
|
||
joined_spans <-st_read(here(data_dir, "spans2.geojson")) %>% st_transform(crs(AllPivots0))
|
||
|
||
pivots_dates <- readRDS(here(harvest_dir, "harvest_data_new")) %>% filter(
|
||
pivot %in% c("1.1", "1.2", "1.3", "1.4", "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") #without 1.6
|
||
)
|
||
AllPivots <- merge(AllPivots0, pivots_dates, by = "pivot_quadrant") %>%
|
||
rename(pivot = pivot.x) #%>% select(-pivot.y)
|
||
|
||
AllPivots_merged <- AllPivots %>%
|
||
group_by(pivot) %>% summarise()
|
||
|
||
AllPivots_merged <- st_transform(AllPivots_merged, crs = proj4string(CI))
|
||
|
||
AllPivots_merged$pivot <- as.factor(AllPivots_merged$pivot)
|
||
AllPivots_merged$pivot <- ordered(AllPivots_merged$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"))
|
||
|
||
pivot_names <- unique(CI_quadrant$pivot)
|
||
|
||
```
|
||
|
||
|
||
```{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){
|
||
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) +
|
||
tm_shape(pivot_shape) +
|
||
tm_borders(lwd = 3) + tm_text("pivot_quadrant", size = 1/2) +
|
||
tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha=0.5)
|
||
}
|
||
|
||
create_CI_diff_map <- function(pivot_raster, pivot_shape, pivot_spans, show_legend = F, legend_is_portrait = F, week_1, week_2, age){
|
||
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) +
|
||
tm_shape(pivot_shape) +
|
||
tm_borders(lwd = 3) + tm_text("pivot_quadrant", size = 1/2) +
|
||
tm_shape(pivot_spans) + tm_borders(lwd = 0.5, alpha=0.5)
|
||
}
|
||
|
||
ci_plot <- function(pivotName){
|
||
# pivotName = "1.1"
|
||
pivotShape <- AllPivots_merged %>% terra::subset(pivot %in% pivotName) %>% st_transform(crs(CI))
|
||
age <- AllPivots %>% dplyr::filter(pivot %in% pivotName) %>% st_drop_geometry() %>% dplyr::select(Age) %>% unique()
|
||
|
||
AllPivots2 <- AllPivots %>% dplyr::filter(pivot %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 <- pivots_dates %>% dplyr::filter(pivot %in% pivotName) %>% ungroup() %>% dplyr::select(planting_date) %>% unique()
|
||
|
||
joined_spans2 <- joined_spans %>% st_transform(crs(pivotShape)) %>% dplyr::filter(pivot %in% pivotName) %>% 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)
|
||
CImap_m1 <- create_CI_map(singlePivot_m1, AllPivots2, joined_spans2, show_legend= F, legend_is_portrait = F, week = week_minus_1, age = age -1)
|
||
CImap <- create_CI_map(singlePivot, AllPivots2, joined_spans2, show_legend= F, legend_is_portrait = F, week = week_minus_1, age = age )
|
||
|
||
|
||
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)
|
||
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)
|
||
|
||
tst <- tmap_arrange(CImap_m2, CImap_m1, CImap,CI_max_abs_last_week, CI_max_abs_three_week, nrow = 1)
|
||
|
||
cat(paste("## Pivot", pivotName, "-", age$Age[1], "weeks after planting/harvest", "\n"))
|
||
# cat("\n")
|
||
# cat('<h2> Pivot', pivotName, '- week', week, '-', age$Age, 'weeks after planting/harvest <h2>')
|
||
# cat(paste("# Pivot",pivots$pivot[i],"\n"))
|
||
print(tst)
|
||
|
||
}
|
||
|
||
|
||
|
||
|
||
cum_ci_plot <- function(pivotName){
|
||
|
||
# pivotName = "1.17"
|
||
data_ci <- CI_quadrant %>% filter(pivot == pivotName)
|
||
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)
|
||
|
||
perfect_pivot_raw <- CI_quadrant %>% group_by(pivot) %>% filter(pivot == "5.1" & season == "2022") %>%
|
||
group_by(DOY) %>% summarise(cumulative_CI = mean(cumulative_CI)) %>% mutate(CI_rate = cumulative_CI/DOY) %>%
|
||
mutate(mean_rolling10 = rollapplyr(CI_rate , width = 10, FUN = mean, partial = TRUE))
|
||
|
||
unique_seasons <- unique(date_preperation_perfect_pivot$season)
|
||
|
||
if(length(unique_seasons) == 3) {
|
||
unique_seasons <- unique_seasons[c(2,3)]
|
||
} else {
|
||
unique_seasons <- unique_seasons
|
||
}
|
||
|
||
perfect_pivot <- perfect_pivot_raw
|
||
|
||
for (s in unique_seasons) {
|
||
season_dates <- seq(from = date_preperation_perfect_pivot$min_date[date_preperation_perfect_pivot$season == s],
|
||
to = date_preperation_perfect_pivot$min_date[date_preperation_perfect_pivot$season == s] + nrow(perfect_pivot_raw) - 1,
|
||
by = "1 day")
|
||
col_name <- as.character(s)
|
||
perfect_pivot <- dplyr::bind_cols(perfect_pivot, tibble(!!col_name := season_dates))
|
||
}
|
||
|
||
perfect_pivot <- perfect_pivot %>%
|
||
pivot_longer(cols = -c("DOY", "cumulative_CI", "CI_rate", "mean_rolling10"),
|
||
names_to = "season", values_to = "Date")
|
||
|
||
g <- ggplot() +
|
||
facet_wrap(~season, scales = "free_x") +
|
||
geom_line(data= data_ci2 %>% filter(season %in% unique_seasons), aes(Date, mean_rolling10, col = Field)) +
|
||
geom_line(data= perfect_pivot, aes(Date , mean_rolling10, col = "Model CI (p5.1 Data 2022, \n date x axis is fictive)"), lty="11",size=1) +
|
||
labs(title = paste("14 day rolling MEAN CI rate - Pivot ", pivotName))+
|
||
# scale_y_continuous(limits=c(0.5,3), breaks = seq(0.5, 3, 0.5))+
|
||
scale_x_date(date_breaks = "1 month", date_labels = "%m-%Y") +
|
||
theme(axis.text.x = element_text(angle = 60, hjust = 1))
|
||
# options(repr.plot.width = 2, repr.plot.height = 2)
|
||
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
|
||
|
||
|
||
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(AllPivots)+ tm_borders( col = "white") +
|
||
tm_text("pivot_quadrant", size = .6, col = "white")
|
||
|
||
```
|
||
\newpage
|
||
|
||
```{r 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(AllPivots)+ tm_borders( col = "black") +
|
||
tm_text("pivot_quadrant", size = .6, col = "black")
|
||
|
||
|
||
|
||
# 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(AllPivots)+ tm_borders( col = "black") +
|
||
# tm_text("pivot_quadrant", size = .6, col = "black")
|
||
|
||
|
||
|
||
# tm_shape(last_week_dif_raster, unit = "m")+
|
||
# tm_raster(breaks = c(-Inf,-50,-25,-5,5,25, Inf), palette = "RdYlGn", midpoint = NA,legend.is.portrait = T) +
|
||
# tm_layout(legend.outside = TRUE,legend.outside.position = "right",legend.show = F, 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(AllPivots)+ tm_borders( col = "black") +
|
||
# tm_text("pivot_quadrant", size = .6, col = "black")
|
||
```
|
||
# Estate fields
|
||
\newpage
|
||
|
||
```{r plots_ci_estate, echo=FALSE, fig.height=3.8, fig.width=10, message=FALSE, warning=FALSE, results='asis'}
|
||
# # pivots <- AllPivots_merged %>% filter(pivot != c("1.1", "1.17"))
|
||
pivots_estate <- AllPivots_merged %>% filter(pivot %in% 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" , "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")
|
||
|
||
walk(pivots_estate$pivot, ~ {
|
||
cat("\n") # Add an empty line for better spacing
|
||
ci_plot(.x)
|
||
cum_ci_plot(.x)
|
||
})
|
||
```
|
||
|
||
|
||
# Coop fields
|
||
\newpage
|
||
|
||
```{r plots_ci_coops, echo=FALSE, fig.height=3.8, fig.width=10, message=FALSE, warning=FALSE, results='asis'}
|
||
pivots_coop <- AllPivots_merged %>% filter(pivot %in% c("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"))
|
||
# pivots_coop <- AllPivots_merged %>% filter(pivot %in% c("2.1", "2.2"))
|
||
|
||
walk(pivots_coop$pivot, ~ {
|
||
cat("\n") # Add an empty line for better spacing
|
||
ci_plot(.x)
|
||
cum_ci_plot(.x)
|
||
})
|
||
|
||
```
|
||
|
||
```{r eval=FALSE, fig.height=10, fig.width=14, include=FALSE}
|
||
CI_all2 <- readRDS(here(cumulative_CI_vals_dir, "All_pivots_Cumulative_CI_whole_pivot_year.rds")) %>%
|
||
mutate(#line = substr(pivot, 1 , 1),
|
||
season = as.factor(season))
|
||
|
||
pivots_dates <- readRDS(here(harvest_dir, "harvest_data_new")) #%>%
|
||
|
||
pvt_age <- pivots_dates %>% ungroup() %>% select(pivot, Age) %>% unique()
|
||
CI_all2 <- left_join(CI_all2, pvt_age , by = "pivot") %>% mutate(month = plyr::round_any((Age/4),2)) %>%
|
||
mutate(month = case_when(month > 16 ~ 18,
|
||
TRUE ~month)) %>%
|
||
group_by(pivot) %>% filter(Age == max(Age)) %>% ungroup()
|
||
|
||
CI_all2$season <- ordered(CI_all2$season, levels = c("Data_2021", "Data_2022"))
|
||
# CI_all2 <- CI_all2 %>% mutate(season_order = as.integer(season))
|
||
latest_model <- CI_all2 %>% group_by(pivot) %>% filter(season =="Data_2022")
|
||
# latest_model <- CI_all2 %>% group_by(pivot) %>% arrange(season, desc(DOY)) %>% slice(1)
|
||
|
||
# CI_all2 <- CI_all %>% group_by(pivot, DOY ) %>% mutate(pivot_cumulative_CI = mean(cumulative_CI))
|
||
# label_data <- CI_all2 %>% group_by(pivot) %>% arrange(season, desc(cumulative_CI)) %>% slice(1) %>% mutate(label = paste(pivot, " - ", round(cumulative_CI)))
|
||
label_data <- latest_model %>% arrange(season, desc(cumulative_CI)) %>% slice(1) %>% mutate(label = paste(pivot, " - ", round(cumulative_CI)))
|
||
|
||
max_day <- label_data %>% group_by(pivot) %>% summarise(max_day = max(DOY))
|
||
|
||
ggplot(data= CI_all2%>% filter(season =="Data_2022"), aes(DOY, cumulative_CI, col = pivot)) +
|
||
facet_wrap(~month) +
|
||
geom_line() +
|
||
# scale_y_continuous(breaks = seq(0, max(label_data$cumulative_CI) + 100, by = 100)) +
|
||
labs(title = "Cumulative CI values over the year per pivot, split per 2-month age (rounded down)", x = "Days after harvest/planting (scale is per 2 weeks)", y = "Cumulative CI value",
|
||
color = "Line") +
|
||
geom_label_repel(data = label_data %>% filter(model %in% latest_model$model), aes(DOY, cumulative_CI, label = label), box.padding = 1,
|
||
# ylim = c(1300, NA),
|
||
max.overlaps = Inf
|
||
# segment.linetype = 4,
|
||
|
||
# segment.curvature = -1e-20,
|
||
# arrow = arrow(length = unit(0.015, "npc"))
|
||
) +
|
||
theme(legend.position="right", legend.text = element_text(size=8), legend.title = element_text(size = 8),
|
||
plot.title = element_text(size=19)) +
|
||
# geom_smooth()+
|
||
guides(fill = guide_legend(nrow=2,byrow=TRUE)) +
|
||
scale_y_continuous(breaks=seq(0,max(label_data$cumulative_CI),100)) +
|
||
scale_x_continuous(breaks=seq(0,max(max_day$max_day),30)) + theme(axis.text.x = element_text(angle = 90)) +
|
||
labs(x = "Weeks")+
|
||
theme(legend.position = "none")
|
||
|
||
```
|
||
# 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")) %>%
|
||
rename( pivot_quadrant = Field)#All_pivots_Cumulative_CI.rds
|
||
ggplot(CI_quadrant %>% filter(pivot %in% "1.11")) +
|
||
geom_line(aes(DOY, cumulative_CI, col = as.factor(season))) +
|
||
facet_wrap(~pivot_quadrant)
|
||
|
||
pivots_dates0 <- readRDS(here(harvest_dir, "harvest_data_new")) %>% ungroup() %>% unique() %>%
|
||
dplyr::select(pivot, pivot_quadrant, Tcha_2021, Tcha_2022 ) %>%
|
||
pivot_longer(cols = c("Tcha_2021", "Tcha_2022"), names_to = "Tcha_Year", values_to = "Tcha") %>%
|
||
filter(Tcha > 50) %>%
|
||
mutate(season = as.integer(str_extract(Tcha_Year, "\\d+")))
|
||
|
||
CI_and_yield <- left_join(CI_quadrant , pivots_dates0, by = c("pivot", "pivot_quadrant", "season")) %>% filter(!is.na(Tcha)) %>%
|
||
group_by(pivot_quadrant, season) %>% slice(which.max(DOY)) %>%
|
||
dplyr::select(pivot, pivot_quadrant, Tcha_Year, Tcha, cumulative_CI, DOY, season) %>%
|
||
mutate(CI_per_day = cumulative_CI/ DOY)
|
||
|
||
ggplot(CI_and_yield) +
|
||
geom_point(aes(Tcha, CI_per_day, col = Tcha_Year ))
|
||
|
||
|
||
set.seed(20)
|
||
CI_and_yield_split <- initial_split(CI_and_yield, prop = 0.75, strata = pivot_quadrant)
|
||
CI_and_yield_test <- training(CI_and_yield_split)
|
||
CI_and_yield_validation <- testing(CI_and_yield_split)
|
||
|
||
|
||
predictors <- c( "cumulative_CI" , "DOY" ,"CI_per_day" )
|
||
response <- "Tcha"
|
||
CI_and_yield_test <- as.data.frame(CI_and_yield_test)
|
||
|
||
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
|
||
)
|
||
|
||
pred_ffs_rf <-
|
||
predict(model_ffs_rf, newdata = CI_and_yield_validation) %>% as.data.frame() %>% rename(predicted_Tcha = ".") %>% mutate(
|
||
pivot_quadrant = CI_and_yield_validation$pivot_quadrant,
|
||
pivot = CI_and_yield_validation$pivot,
|
||
Age_days = CI_and_yield_validation$DOY,
|
||
total_CI = round(CI_and_yield_validation$cumulative_CI, 0),
|
||
predicted_Tcha = round(predicted_Tcha, 0),
|
||
season = CI_and_yield_validation$season
|
||
) %>% dplyr::select(pivot , pivot_quadrant, Age_days, total_CI, predicted_Tcha, season) %>%
|
||
left_join(., CI_and_yield_validation, by = c("pivot", "pivot_quadrant", "season")) %>%
|
||
filter(Age_days > 250)
|
||
|
||
|
||
|
||
prediction_2023 <- CI_quadrant %>% filter(season == "2023") %>% group_by(pivot_quadrant) %>% slice(which.max(DOY))%>%
|
||
mutate(CI_per_day = cumulative_CI/ DOY)
|
||
|
||
pred_rf_2023 <- predict(model_ffs_rf, newdata=prediction_2023) %>%
|
||
as.data.frame() %>% rename(predicted_Tcha_2023 = ".") %>% mutate(pivot_quadrant = prediction_2023$pivot_quadrant,
|
||
pivot = prediction_2023$pivot,
|
||
Age_days = prediction_2023$DOY,
|
||
total_CI = round(prediction_2023$cumulative_CI,0),
|
||
predicted_Tcha_2023 = round(predicted_Tcha_2023, 0)) %>%
|
||
filter(Age_days > 300) %>%
|
||
dplyr::select(pivot ,pivot_quadrant, Age_days, total_CI, predicted_Tcha_2023)%>%
|
||
mutate(CI_per_day = round(total_CI/ Age_days, 1))
|
||
|
||
|
||
|
||
```
|
||
|
||
```{r echo=FALSE}
|
||
ggplot(pred_ffs_rf, aes(y = predicted_Tcha , x = Tcha , col = pivot )) +
|
||
geom_point() +geom_abline() +
|
||
scale_x_continuous(limits = c(50, 160))+
|
||
scale_y_continuous(limits = c(50, 160)) +
|
||
labs(title = "Model trained and tested on historical results - RF")
|
||
|
||
ggplot(pred_rf_2023, aes(total_CI , predicted_Tcha_2023 , col = pivot )) +
|
||
geom_point() + labs(title = "2023 data (still to be harvested) - Fields over 300 days old")
|
||
|
||
knitr::kable(pred_rf_2023)
|
||
```
|