--- # title: paste0("CI report week ", week, " - all pivots from ", last_tuesday, " to ", today) params: ref: word-styles-reference-03.docx output: word_document: reference_docx: !expr file.path("word-styles-reference-03.docx") # toc: true editor_options: chunk_output_type: console --- ```{r setup, include=FALSE} # 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} library(here) library(sf) library(tidyverse) library(tmap) library(lubridate) library(exactextractr) library(zoo) library(raster) library(rsample) library(caret) library(CAST) ``` ```{r directories, message=FALSE, warning=FALSE, include=FALSE} laravel_storage_dir <- here("laravel_app/storage/app") 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(data_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-10-05" week <- week(today) #today = "2022-08-18" #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") ``` --- title: `r title_var` --- 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=FALSE, warning=FALSE, include=FALSE} # get latest CI index week_minus_1 <- as.numeric(week) -1 week_minus_2 <- as.numeric(week) -2 week_minus_3 <- as.numeric(week) -3 today_minus_1 <- as.character(ymd(today) - 7) today_minus_2 <- as.character(ymd(today) - 14) today_minus_3 <- as.character(ymd(today) - 21) # 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, "_2023.tif"))) %>% subset("CI") CI_m1 <- brick(here(weekly_CI_mosaic, paste0("week_",week_minus_1, "_2023.tif"))) %>% subset("CI") #%>% subset("CI") CI_m2 <- brick(here(weekly_CI_mosaic, paste0("week_",week_minus_2, "_2023.tif"))) %>% subset(6) #%>% subset("CI") CI_m3 <- brick(here(weekly_CI_mosaic, paste0("week_",week_minus_3, "_2023.tif"))) %>% subset(6) #%>% subset("CI") last_week_dif_raster <- ((CI - CI_m1) / CI_m1) * 100 last_week_dif_raster_abs <- (CI - CI_m1) two_week_dif_raster_abs <- (CI - CI_m2) 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)) } ci_plot <- function(pivotName){ # pivotName = "2.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_two_week <- two_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_m1, 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_two_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_2, age = age) tst <- tmap_arrange(CImap_m2, CImap_m1, CImap,CI_max_abs_last_week, CI_max_abs_two_week, nrow = 1) cat('

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

') print(tst) } 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) } 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(pivot_quadrant) %>% 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 = pivot_quadrant)) + 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 eval=FALSE, fig.height=7.2, fig.width=10, message=FALSE, warning=FALSE, include=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 = 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") 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") ``` ```{r plots_ci, echo=FALSE, fig.height=3.7, fig.width=10, message=FALSE, warning=FALSE, results='asis'} # ci_plot("1.17") # cum_ci_plot("1.17") # x = 1 # for(j in x){ # coops <- Chemba_pivot_owners %>% filter(`OWNER update 18/6/2022` %in% c("chapo", "Lambane", "Canhinbe" )) pivots <- AllPivots_merged %>% filter(pivot != "1.17") #%>% filter(pivot %in% c( "2.1", "2.2", "2.3", "2.4", "3.1", "3.2", "3.3", "4.4", "4.6" , "4.3", "4.5", "4.2", "4.1", "5.1", "5.2", "5.3", "5.4", "7.1", "7.2", "7.3" , "7.4", "7.5", "7.6" )) for(i in pivots$pivot) { ci_plot(i) cum_ci_plot(i) } # lapply(pivots, function(pivot) { # ci_plot(pivot) # cum_ci_plot(pivot) # }) #cat("\\newpage") # } ``` ```{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") ``` The below table shows estimates of the biomass if you would harvest them now. ```{r eval=FALSE, 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 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) CI_and_yield <- left_join(CI_quadrant , pivots_dates0, by = c("pivot", "pivot_quadrant")) %>% filter(!is.na(Tcha)) %>% group_by(pivot_quadrant) %>% slice(which.max(DOY)) %>% dplyr::select(pivot, pivot_quadrant, Tcha_Year, Tcha, cumulative_CI, DOY) %>% mutate(CI_per_day = cumulative_CI/ DOY) set.seed(20) CI_and_yield_split <- initial_split(CI_and_yield, prop = 0.75, weight = 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) ) %>% dplyr::select(pivot , pivot_quadrant, Age_days, total_CI, predicted_Tcha) %>% left_join(., CI_and_yield_validation, by = c("pivot", "pivot_quadrant")) %>% filter(Age_days > 250) 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") prediction_2023 <- CI_quadrant %>% filter(season == "Data_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)) 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) ``` ```{r eval=FALSE, include=FALSE} model_CI <-lm( formula = cumulative_CI ~ DOY , data = CI_and_yield_test ) pivot_ = "4.4" df4 = data.frame(pivot_, 365, NA) names(df4)=c("pivot", "DOY", "cumulative_CI") a <- CI_all %>% filter(season == "Data_2022", pivot == pivot_) %>% ungroup() %>% select(pivot, DOY, cumulative_CI) %>% complete(DOY = seq.int(max(DOY), 365, 1), pivot = pivot_) %>% arrange(DOY) # complete(DOY = seq.int(max(DOY), 365, 1)) # rbind(.,df4) b <- predict(model_CI, a) %>% as.data.frame() %>% slice(which.max(.)) %>% rename(cumulative_CI = ".") %>% mutate(DOY = 365) pred_CI_2022 <- predict(model, newdata=b ) %>% as.data.frame() %>% rename(predicted_Tcha_365 = ".") %>% mutate(pivot = df4$pivot, predicted_Tcha_365 = round(predicted_Tcha_365, 0), Age_days = df4$DOY) pred_CI_2022 ```