144 lines
4.3 KiB
R
144 lines
4.3 KiB
R
|
|
library(here)
|
|
library(sf)
|
|
library(terra)
|
|
library(tidyverse)
|
|
library(lubridate)
|
|
library(exactextractr)
|
|
# library(readxl)
|
|
|
|
|
|
# Vang alle command line argumenten op
|
|
args <- commandArgs(trailingOnly = TRUE)
|
|
|
|
|
|
# Converteer het tweede argument naar een string waarde
|
|
project_dir <- as.character(args[1])
|
|
|
|
# Controleer of data_dir een geldige waarde is
|
|
if (!is.character(project_dir)) {
|
|
project_dir <- "chemba"
|
|
}
|
|
|
|
|
|
laravel_storage_dir <- here("laravel_app/storage/app", project_dir)
|
|
#preparing directories
|
|
planet_tif_folder <- here(laravel_storage_dir, "merged_tif")
|
|
merged_final <- here(laravel_storage_dir, "merged_final_tif")
|
|
|
|
new_project_question = FALSE
|
|
|
|
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")
|
|
|
|
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"))
|
|
# source(here("r_app/ci_extraction_utils.R"))
|
|
# source(here("r_app/mosaic_creation_utils.R"))
|
|
|
|
source("parameters_project.R")
|
|
source("ci_extraction_utils.R")
|
|
|
|
pivot_stats2 <- readRDS(here(cumulative_CI_vals_dir,"combined_CI_data.rds")) %>%
|
|
ungroup() %>%
|
|
group_by(field, sub_field) %>%
|
|
summarise(across(everything(), ~ first(na.omit(.))), .groups = "drop")
|
|
|
|
#%>% drop_na(pivot_quadrant)
|
|
|
|
# gather data into long format for easier calculation and visualisation
|
|
pivot_stats_long <- pivot_stats2 %>%
|
|
tidyr::gather("Date", value, -field, -sub_field ) %>%
|
|
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,
|
|
# sub_field = field) %>%
|
|
unique()
|
|
|
|
years <- harvesting_data %>%
|
|
filter(!is.na(season_start)) %>%
|
|
distinct(year) %>%
|
|
pull(year)
|
|
|
|
extract_CI_data <- function(field_names, harvesting_data, field_CI_data, season) {
|
|
# Filter harvesting data for the given season and field names
|
|
filtered_harvesting_data <- harvesting_data %>%
|
|
filter(year == season, sub_field %in% field_names)
|
|
|
|
# Filter field CI data for the given field names
|
|
filtered_field_CI_data <- field_CI_data %>%
|
|
filter(sub_field %in% field_names)
|
|
|
|
# Return an empty data frame if no CI data is found
|
|
if (nrow(filtered_field_CI_data) == 0) {
|
|
return(data.frame())
|
|
}
|
|
|
|
# Create a linear interpolation function for the CI data
|
|
ApproxFun <- approxfun(x = filtered_field_CI_data$Date, y = filtered_field_CI_data$value)
|
|
Dates <- seq.Date(ymd(min(filtered_field_CI_data$Date)), ymd(max(filtered_field_CI_data$Date)), by = 1)
|
|
LinearFit <- ApproxFun(Dates)
|
|
|
|
# Combine interpolated data with the original CI data
|
|
CI <- data.frame(Date = Dates, FitData = LinearFit) %>%
|
|
left_join(., filtered_field_CI_data, by = "Date") %>%
|
|
filter(Date > filtered_harvesting_data$season_start & Date < filtered_harvesting_data$season_end)
|
|
|
|
# If CI is empty after filtering, return an empty dataframe
|
|
if (nrow(CI) == 0) {
|
|
return(data.frame())
|
|
}
|
|
|
|
# Add additional columns if data exists
|
|
CI <- CI %>%
|
|
mutate(DOY = seq(1, n(), 1),
|
|
model = paste0("Data", season, " : ", field_names),
|
|
season = season,
|
|
subField = field_names)
|
|
|
|
return(CI)
|
|
}
|
|
|
|
|
|
|
|
CI_all <- map_df(years, function(yr) {
|
|
# yr = 2021
|
|
message(yr)
|
|
|
|
# Get the fields harvested in this year
|
|
sub_fields <- harvesting_data %>%
|
|
filter(year == yr) %>%
|
|
filter(!is.na(season_start)) %>%
|
|
pull(sub_field)
|
|
|
|
# Filter sub_fields to only include those with value data in pivot_stats_long
|
|
valid_sub_fields <- sub_fields %>%
|
|
keep(~ any(pivot_stats_long$sub_field == .x))
|
|
|
|
# Extract data for each valid field
|
|
map(valid_sub_fields, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = yr)) %>%
|
|
list_rbind()
|
|
})
|
|
|
|
CI_all <- CI_all %>% group_by(model) %>% mutate(CI_per_day = FitData - lag(FitData),
|
|
cumulative_CI = cumsum(FitData))
|
|
|
|
|
|
message('CI_all cumulative')
|
|
|
|
head(CI_all)
|
|
message('show head')
|
|
|
|
saveRDS(CI_all, here(cumulative_CI_vals_dir,"All_pivots_Cumulative_CI_quadrant_year_v2.rds"))
|
|
message('rds saved')
|
|
|