adjusted scripts to be more generalised

This commit is contained in:
Timon 2024-03-13 10:27:10 +01:00
parent 238520ff89
commit 66d6f80532
7 changed files with 1412 additions and 247 deletions

1
r_app/.gitignore vendored
View file

@ -16,4 +16,5 @@ renv
.Rprofile
# Ignore OSX files
.DS_Store
CI_report_dashboard_planet_files

View file

@ -11,7 +11,6 @@ library(lubridate)
library(exactextractr)
library(CIprep)
source(here("r_app", "parameters_project.R"))
# Vang alle command line argumenten op
args <- commandArgs(trailingOnly = TRUE)
@ -57,6 +56,8 @@ weekly_CI_mosaic <- here(laravel_storage_dir, "weekly_mosaic")
daily_vrt <- here(data_dir, "vrt")
harvest_dir <- here(data_dir, "HarvestData")
source(here("r_app", "parameters_project.R"))
dir.create(here(laravel_storage_dir))
dir.create(here(data_dir))
dir.create(here(extracted_CI_dir))
@ -134,6 +135,12 @@ create_mask_and_crop <- function(file, field_boundaries) {
vrt_list <- list()
# for (file in raster_files) {
# v_crop <- create_mask_and_crop(file, field_boundaries)
# message(file, " processed")
# gc()
# }
for (file in filtered_files) {
v_crop <- create_mask_and_crop(file, field_boundaries)
emtpy_or_full <- global(v_crop, "notNA")
@ -271,14 +278,14 @@ if (new_project_question == TRUE) {
pivot_stats <- extracted_values %>%
map(readRDS) %>% list_rbind() %>%
group_by(pivot_quadrant) %>%
group_by(subField) %>%
summarise(across(everything(), ~ first(na.omit(.))))
saveRDS(pivot_stats, here(cumulative_CI_vals_dir,"combined_CI_data.rds")) #used to save the rest of the data into one file
pivot_stats2 <- pivot_stats
print("All CI values extracted from allhistoric images")
print("All CI values extracted from all historic images")
} else {
print("combined_CI_data.rds exists, adding the latest image data to the table.")
@ -309,15 +316,15 @@ if (new_project_question == TRUE) {
# gather data into long format for easier calculation and visualisation
pivot_stats_long <- pivot_stats2 %>%
tidyr::gather("Date", value, -pivot_quadrant, -pivot ) %>%
tidyr::gather("Date", value, -Field, -subField, -sub_area ) %>%
mutate(Date = right(Date, 8),
Date = lubridate::ymd(Date)
) %>%
drop_na(c("value","Date")) %>%
mutate(value = as.numeric(value))%>%
filter_all(all_vars(!is.infinite(.))) %>%
rename(Field = pivot_quadrant,
subField = Field) %>%
# rename(Field = pivot_quadrant,
# subField = Field) %>%
unique()

File diff suppressed because one or more lines are too long

View file

@ -1,27 +1,28 @@
---
# 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"
ref: "word-styles-reference-var1.docx"
output_file: CI_report.docx
report_date: "2023-12-12"
data_dir: "chemba"
output:
html_document:
toc: yes
df_print: paged
word_document:
reference_docx: !expr file.path("word-styles-reference-var1.docx")
toc: true
reference_docx: "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
# Activeer de renv omgeving
renv::activate()
# 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
@ -38,7 +39,7 @@ library(lubridate)
library(exactextractr)
library(zoo)
library(raster)
library(terra)
library(rsample)
library(caret)
library(randomForest)
@ -49,15 +50,16 @@ library(CAST)
```
```{r directories, message=FALSE, warning=FALSE, include=FALSE}
laravel_storage_dir <- here("laravel_app/storage/app/",params$data_dir)
project_dir <- "chemba"
# laravel_storage_dir <- here("laravel_app/storage/app/",params$data_dir)
laravel_storage_dir <- here("laravel_app/storage/app",project_dir)
data_dir_project <- here(laravel_storage_dir, "Data")
# message('DATA_DIR',data_dir_project)
extracted_CI_dir <- here(data_dir_project, "extracted_ci")
data_dir <- here(laravel_storage_dir, "Data")
# message('DATA_DIR',data_dir)
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_project, "HarvestData")
harvest_dir <- here(data_dir, "HarvestData")
weekly_CI_mosaic <- here(laravel_storage_dir, "weekly_mosaic")
@ -69,6 +71,7 @@ source(here("r_app", "parameters_project.R"))
```{r week, message=FALSE, warning=FALSE, include=FALSE}
# week <- 5
#today = "2023-12-12"
report_date <- Sys.Date()
today <- as.character(report_date)
week <- week(today)
@ -135,14 +138,14 @@ AllPivots0 <- field_boundaries_sf
AllPivots <- merge(AllPivots0, harvesting_data, by = c("Field", "subField")) #%>%
# rename(pivot = pivot.x) #%>% select(-pivot.y)
head(AllPivots)
head(AllPivots)
AllPivots_merged <- AllPivots %>%
group_by(Field) %>% summarise()
group_by(Field) %>% summarise(sub_area = first(sub_area))
AllPivots_merged <- st_transform(AllPivots_merged, crs = proj4string(CI))
pivot_names <- unique(CI_quadrant$pivot)
pivot_names <- unique(CI_quadrant$Field)
```
@ -228,7 +231,7 @@ create_CI_diff_map <- function(pivot_raster, pivot_shape, pivot_spans, show_lege
}
ci_plot <- function(pivotName){
# pivotName = "1.1"
# pivotName = "AG1D06P"
pivotShape <- AllPivots_merged %>% terra::subset(Field %in% pivotName) %>% st_transform(crs(CI))
age <- AllPivots %>% dplyr::filter(Field %in% pivotName) %>% st_drop_geometry() %>% dplyr::select(Age) %>% unique()
@ -243,9 +246,9 @@ ci_plot <- function(pivotName){
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(Field %in% pivotName) %>% ungroup() %>% dplyr::select(planting_date) %>% unique()
planting_date <- harvesting_data %>% dplyr::filter(Field %in% pivotName) %>% ungroup() %>% dplyr::select(Season_start) %>% unique()
joined_spans2 <- joined_spans %>% st_transform(crs(pivotShape)) %>% dplyr::filter(Field %in% pivotName) %>% unique() %>% st_crop(., pivotShape)
joined_spans2 <- joined_spans %>% st_transform(crs(pivotShape)) %>% dplyr::filter(Field %in% pivotName) #%>% unique() %>% 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)
@ -270,7 +273,7 @@ ci_plot <- function(pivotName){
cum_ci_plot <- function(pivotName){
# pivotName = "1.17"
# pivotName = "1.1"
data_ci <- CI_quadrant %>% filter(Field == pivotName)
data_ci2 <- data_ci %>% mutate(CI_rate = cumulative_CI/DOY,
week = week(Date))%>% group_by(Field) %>%
@ -371,10 +374,9 @@ tm_shape(CI, unit = "m")+
```
# Estate fields
\newpage
```{r plots_ci_estate, echo=FALSE, fig.height=3.8, fig.width=10, message=FALSE, warning=FALSE, results='asis'}
```{r plots_ci_estate, eval=FALSE, fig.height=3.8, fig.width=10, message=FALSE, warning=FALSE, include=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")
@ -397,7 +399,7 @@ pivots_grouped <- AllPivots_merged %>%
# Iterate over each subgroup
for (subgroup in unique(pivots_grouped$sub_area)) {
cat("\n## Subgroup:", subgroup, "\n") # Add a title for the subgroup
cat("\n# Subgroup: ", subgroup, "\n") # Add a title for the subgroup
subset_data <- filter(pivots_grouped, sub_area == subgroup)
walk(subset_data$Field, ~ {
cat("\n") # Add an empty line for better spacing
@ -457,15 +459,15 @@ ggplot(data= CI_all2%>% filter(season =="Data_2022"), aes(DOY, cumulative_CI, co
```
# Yield prediction
The below table shows estimates of the biomass if you would harvest them now.
```{r message=FALSE, warning=FALSE, include=FALSE}
```{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
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 ) %>%
pivots_dates0 <- pivots_dates0 %>% ungroup() %>% unique() %>%
dplyr::select(Field, subField, 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+")))
@ -537,7 +539,7 @@ pred_rf_2023 <- predict(model_ffs_rf, newdata=prediction_2023) %>%
```
```{r yield_plaatjes, echo=FALSE}
```{r yield_plaatjes, eval=FALSE, include=FALSE}
ggplot(pred_ffs_rf, aes(y = predicted_Tcha , x = Tcha , col = pivot )) +
geom_point() +geom_abline() +
scale_x_continuous(limits = c(50, 160))+

Binary file not shown.

After

Width:  |  Height:  |  Size: 136 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 161 KiB

View file

@ -6,8 +6,7 @@ if(project_dir == "chemba"){
field_boundaries_sf <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant)%>%
mutate(sub_area = case_when(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") ~ "estate",
mutate(sub_area = case_when(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") ~ "estate",
TRUE ~ "Cooperative"))
names(field_boundaries_sf) <- c("Field", "subField", "geometry", "sub_area")
@ -16,7 +15,7 @@ if(project_dir == "chemba"){
# field_boundaries_merged <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% group_by(pivot) %>% summarise() %>% vect()
joined_spans <-st_read(here(data_dir_project, "span.geojson")) %>% st_transform(crs(field_boundaries_sf))
joined_spans <-st_read(here(data_dir, "span.geojson")) %>% st_transform(crs(field_boundaries_sf))
names(joined_spans) <- c("Field", "area", "radius", "spans", "span", "geometry")
@ -27,14 +26,14 @@ if(project_dir == "chemba"){
)
harvesting_data <- pivots_dates0 %>%
dplyr::select(c("pivot_quadrant", "season_start_2021", "season_end_2021", "season_start_2022",
dplyr::select(c("pivot_quadrant", "pivot", "season_start_2021", "season_end_2021", "season_start_2022",
"season_end_2022", "season_start_2023", "season_end_2023", "season_start_2024", "season_end_2024", "Age")) %>%
pivot_longer(cols = starts_with("season"), names_to = "Year", values_to = "value") %>%
separate(Year, into = c("name", "Year"), sep = "(?<=season_start|season_end)\\_", remove = FALSE) %>%
mutate(name = str_to_title(name)) %>%
pivot_wider(names_from = name, values_from = value) %>%
rename(subField = pivot_quadrant) %>%
mutate(Field = substr(subField, 1, 3))
rename(Field = pivot,
subField = pivot_quadrant)
@ -46,11 +45,12 @@ if(project_dir == "chemba"){
} else if (project_dir == "xinavane"){
library(readxl)
message("Yield data for Xinavane")
field_boundaries_sf <- st_read(here(data_dir, "Xinavane_demo.geojson")) %>% dplyr::select(-Pivot)
names(field_boundaries_sf) <- c("Field", "subField", "geometry")
field_boundaries <- st_read(here(data_dir, "Xinavane_demo.geojson")) %>% vect()
field_boundaries_sf <- st_read(here(data_dir, "Xinavane_demo.geojson"))
field_boundaries <- field_boundaries_sf %>% vect()
joined_spans <- field_boundaries
joined_spans <- field_boundaries_sf %>% dplyr::select(Field)
pivots_dates0 <- read_excel(here(harvest_dir, "Yield data.xlsx"),
skip = 3,
@ -65,11 +65,11 @@ if(project_dir == "chemba"){
tcha = 6,
tchy = 7,
Season_end = 8,
age = 9,
Age = 9,
ratoon = 10
) %>%
mutate(Season_end = ymd(Season_end),
Season_start = as.Date(Season_end - (duration(months = age))),
Season_start = as.Date(Season_end - (duration(months = Age))),
subField = Field) #don't forget to add 2022 as a year for the 'curent' season
pivots_dates0 <- pivots_dates0 %>%
@ -84,7 +84,7 @@ if(project_dir == "chemba"){
mutate(Season_start = Season_start + years(6))
harvesting_data <- pivots_dates0 %>% select(c("Field","subField", "Year", "Season_start","Season_end", "Age" ))
harvesting_data <- pivots_dates0 %>% dplyr::select(c("Field","subField", "Year", "Season_start","Season_end", "Age" , "sub_area", "tcha"))
} else {