diff --git a/r_app/CI_report_dashboard_planet.Rmd b/r_app/CI_report_dashboard_planet.Rmd index fe9f597..88f4bea 100644 --- a/r_app/CI_report_dashboard_planet.Rmd +++ b/r_app/CI_report_dashboard_planet.Rmd @@ -1,12 +1,12 @@ --- # title: paste0("CI report week ", week, " - all pivots from ", last_tuesday, " to ", today) params: - ref: word-styles-reference-03.docx + 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-03.docx") + reference_docx: !expr file.path("word-styles-reference-var1.docx") # toc: true editor_options: chunk_output_type: console @@ -28,6 +28,7 @@ renv::activate() ``` ```{r libraries, message=FALSE, warning=FALSE, include=FALSE} +knitr::opts_chunk$set(warning = FALSE, message = FALSE) library(here) library(sf) library(tidyverse) @@ -39,6 +40,7 @@ library(raster) library(rsample) library(caret) +library(randomForest) library(CAST) ``` @@ -67,24 +69,9 @@ 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") -``` +#title_var <- paste0("CI dashboard week ", week, " - all pivots dashboard using 3x3 meter resolution") +subtitle_var <- paste("Report generated on", today) ---- -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=TRUE, warning=TRUE, include=FALSE} -# get latest CI index today_minus_1 <- as.character(ymd(today) - 7) today_minus_2 <- as.character(ymd(today) - 14) today_minus_3 <- as.character(ymd(today) - 21) @@ -98,6 +85,23 @@ 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) @@ -108,9 +112,9 @@ CI_m1 <- brick(here(weekly_CI_mosaic, paste0("week_",week_minus_1, "_", year_2, 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 <- ((CI - CI_m1) / CI_m1) * 100 last_week_dif_raster_abs <- (CI - CI_m1) -two_week_dif_raster_abs <- (CI - CI_m2) +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)) @@ -196,45 +200,6 @@ subchunkify <- function(g, fig_height=7, fig_width=5) { 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, 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) + @@ -255,6 +220,46 @@ create_CI_diff_map <- function(pivot_raster, pivot_shape, pivot_spans, show_lege 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('

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

') + # cat(paste("# Pivot",pivots$pivot[i],"\n")) + print(tst) + + } + + + cum_ci_plot <- function(pivotName){ @@ -327,7 +332,7 @@ tm_shape(RGB_raster, unit = "m") + tm_rgb(r=1, g=2, b=3, max.value = 255) + ``` \newpage -```{r eval=FALSE, fig.height=7.2, fig.width=10, message=FALSE, warning=FALSE, include=FALSE} +```{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)")+ @@ -339,50 +344,56 @@ tm_shape(CI, unit = "m")+ -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_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_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 - 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_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) +}) ``` -```{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") +# Coop fields +\newpage -#%>% 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" )) +```{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")) -for(i in pivots$pivot) { - ci_plot(i) - cum_ci_plot(i) - } +walk(pivots_coop$pivot, ~ { + cat("\n") # Add an empty line for better spacing + ci_plot(.x) + cum_ci_plot(.x) +}) -# lapply(pivots, function(pivot) { -# ci_plot(pivot) -# cum_ci_plot(pivot) -# }) - - #cat("\\newpage") -# } ``` ```{r eval=FALSE, fig.height=10, fig.width=14, include=FALSE} @@ -433,22 +444,30 @@ ggplot(data= CI_all2%>% filter(season =="Data_2022"), aes(DOY, cumulative_CI, co theme(legend.position = "none") ``` - +# Yield prediction The below table shows estimates of the biomass if you would harvest them now. -```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE} +```{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) + 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")) %>% filter(!is.na(Tcha)) %>% - group_by(pivot_quadrant) %>% slice(which.max(DOY)) %>% - dplyr::select(pivot, pivot_quadrant, Tcha_Year, Tcha, cumulative_CI, DOY) %>% +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) @@ -482,18 +501,15 @@ pred_ffs_rf <- 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")) %>% + 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) -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))%>% + +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) %>% @@ -506,32 +522,19 @@ pred_rf_2023 <- predict(model_ffs_rf, newdata=prediction_2023) %>% 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) - ``` - -```{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 -``` - diff --git a/r_app/Rplots.pdf b/r_app/Rplots.pdf index 2a4ceb7..0755b82 100644 Binary files a/r_app/Rplots.pdf and b/r_app/Rplots.pdf differ diff --git a/runpython.sh b/runpython.sh index 0a0a832..3c395b6 100755 --- a/runpython.sh +++ b/runpython.sh @@ -3,6 +3,8 @@ date=$(date +%Y-%m-%d) # Standaardwaarde voor days days=1 +project_dir="chemba" + # Loop door alle argumenten while [ "$#" -gt 0 ]; do case "$1" in @@ -12,6 +14,9 @@ while [ "$#" -gt 0 ]; do --date=*) date="${1#*=}" ;; + --project_dir=*) + project_dir="${1#*=}" + ;; *) echo "Onbekende optie: $1" exit 1 @@ -30,6 +35,7 @@ source "$script_dir/python_app/myenv/bin/activate" export DAYS=$days export DATE=$date +export ProjectDir=$project_dir # Hier kan je verdere stappen toevoegen, zoals het uitvoeren van je Python-script of Jupyter Notebook jupyter nbconvert --execute --to script --stdout "$script_dir/python_app/Chemba_download.ipynb"