revised report

This commit is contained in:
Martin Folkerts 2024-01-16 14:43:22 +01:00
parent bce3604ae5
commit d3c2c18918
3 changed files with 144 additions and 135 deletions

View file

@ -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('<h1> Pivot', pivotName, '- week', week, '-', age$Age, 'weeks after planting/harvest <h1>')
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('<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){
@ -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
```

Binary file not shown.

View file

@ -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"