wip
This commit is contained in:
Martin Folkerts
parent
915284ba04
commit
1f194c5670
29
interpolate_growth_model.sh
Executable file
29
interpolate_growth_model.sh
Executable file
|
|
@ -0,0 +1,29 @@
|
|||
#!/bin/bash
|
||||
|
||||
project_dir="chemba"
|
||||
|
||||
# Parse command line arguments
|
||||
for arg in "$@"; do
|
||||
case $arg in
|
||||
--project_dir=*)
|
||||
project_dir="${arg#*=}"
|
||||
;;
|
||||
*)
|
||||
echo "Unknown option: $arg"
|
||||
exit 1
|
||||
;;
|
||||
esac
|
||||
shift
|
||||
done
|
||||
|
||||
|
||||
# Check if required arguments are set
|
||||
if [ -z "$project_dir" ]; then
|
||||
echo "Missing argument project_dir. Use: interpolate_growth_model.sh --project_dir=chemba"
|
||||
exit 1
|
||||
fi
|
||||
|
||||
echo interpolate_growth_model.R $project_dir
|
||||
|
||||
cd ../r_app
|
||||
Rscript interpolate_growth_model.R $project_dir
|
||||
|
|
@ -56,6 +56,7 @@ harvest_dir <- here(data_dir, "HarvestData")
|
|||
|
||||
source("parameters_project.R")
|
||||
source("ci_extraction_utils.R")
|
||||
source("mosaic_creation_utils.R")
|
||||
|
||||
dir.create(here(laravel_storage_dir))
|
||||
dir.create(here(data_dir))
|
||||
|
|
@ -86,7 +87,60 @@ print(dates)
|
|||
# flatten_chr()
|
||||
# head(filtered_files)
|
||||
|
||||
raster_files <- list.files(planet_tif_folder,full.names = T, pattern = ".tif")
|
||||
|
||||
filtered_files <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>%
|
||||
compact() %>%
|
||||
flatten_chr()
|
||||
|
||||
# filtered_files <- raster_files #for first CI extraction
|
||||
|
||||
# create_mask_and_crop <- function(file, field_boundaries) {
|
||||
# message("starting ", file)
|
||||
# CI <- rast(file)
|
||||
# # names(CI) <- c("green","nir")
|
||||
# message("raster loaded")
|
||||
#
|
||||
# CI <- CI[[2]]/CI[[1]]-1
|
||||
# # CI <- CI$nir/CI$green-1
|
||||
# message("CI calculated")
|
||||
# CI <- terra::crop(CI, field_boundaries, mask = TRUE) #%>% CI_func()
|
||||
#
|
||||
# new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif"))
|
||||
# writeRaster(CI, new_file, overwrite = TRUE)
|
||||
#
|
||||
# vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
|
||||
# terra::vrt(new_file, vrt_file, overwrite = TRUE)
|
||||
#
|
||||
# return(CI)
|
||||
# }
|
||||
|
||||
|
||||
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")
|
||||
|
||||
vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
|
||||
|
||||
if(emtpy_or_full[1,] > 10000){
|
||||
vrt_list[vrt_file] <- vrt_file
|
||||
|
||||
}else{
|
||||
file.remove(vrt_file)
|
||||
|
||||
}
|
||||
|
||||
message(file, " processed")
|
||||
gc()
|
||||
}
|
||||
|
||||
# Extracting CI
|
||||
|
||||
|
|
@ -129,7 +183,6 @@ print("combined_CI_data.rds exists, adding the latest image data to the table.")
|
|||
filtered_files <- map(dates$days_filter, ~ raster_files_NEW[grepl(pattern = .x, x = raster_files_NEW)]) %>%
|
||||
compact() %>%
|
||||
flatten_chr()
|
||||
|
||||
walk(filtered_files, extract_rasters_daily, field_geojson= field_boundaries, quadrants = TRUE, daily_CI_vals_dir)
|
||||
|
||||
extracted_values <- list.files(daily_CI_vals_dir, full.names = TRUE)
|
||||
|
|
@ -140,103 +193,11 @@ extracted_values <- map(dates$days_filter, ~ extracted_values[grepl(pattern = .x
|
|||
|
||||
pivot_stats <- extracted_values %>%
|
||||
map(readRDS) %>% list_rbind() %>%
|
||||
group_by(sub_field) %>%
|
||||
summarise(across(everything(), ~ first(na.omit(.))))
|
||||
group_by(sub_field)
|
||||
|
||||
|
||||
combined_CI_data <- readRDS(here(cumulative_CI_vals_dir,"combined_CI_data.rds")) #%>% drop_na(pivot_quadrant)
|
||||
pivot_stats2 <- bind_rows(pivot_stats, combined_CI_data)
|
||||
# pivot_stats2 <- combined_CI_data
|
||||
print("All CI values extracted from latest 7 images.")
|
||||
saveRDS(combined_CI_data, here(cumulative_CI_vals_dir,"combined_CI_data.rds")) #used to save the rest of the data into one file
|
||||
|
||||
|
||||
# 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()
|
||||
|
||||
|
||||
# #2021
|
||||
# pivot_select_model_Data_2021 <- harvesting_data %>% filter(Year == 2021) %>% pull(field)
|
||||
#
|
||||
# pivot_select_model_Data_2022 <- harvesting_data %>% filter(Year == 2022) %>% pull(field)
|
||||
|
||||
# pivot_select_model_Data_2023 <- harvesting_data %>% filter(year == 2023) %>% filter(!is.na(season_start)) %>% pull(sub_field)
|
||||
|
||||
pivot_select_model_Data_2024 <- harvesting_data %>% filter(year == 2024)%>% filter(!is.na(season_start)) %>% pull(sub_field)
|
||||
|
||||
# pivots_dates_Data_2022 <- pivots_dates0 %>% filter(!is.na(season_end_2022))
|
||||
# pivot_select_model_Data_2022 <- unique(pivots_dates_Data_2022$pivot_quadrant )
|
||||
#
|
||||
# pivots_dates_Data_2023 <- pivots_dates0 %>% filter(!is.na(season_start_2023))
|
||||
# pivot_select_model_Data_2023 <- unique(pivots_dates_Data_2023$pivot_quadrant)
|
||||
#
|
||||
# pivots_dates_Data_2024 <- pivots_dates0 %>% filter(!is.na(season_start_2024))
|
||||
# pivot_select_model_Data_2024 <- unique(pivots_dates_Data_2024$pivot_quadrant)
|
||||
|
||||
extract_CI_data <- function(field_names, harvesting_data, field_CI_data, season) { # nolint: line_length_linter.
|
||||
# field_names = "AG1D06P"
|
||||
# field_names = "1.13A"
|
||||
# harvesting_data = harvesting_data
|
||||
# field_CI_data = pivot_stats_long
|
||||
# season= 2023
|
||||
|
||||
filtered_harvesting_data <- harvesting_data %>%
|
||||
na.omit() %>%
|
||||
filter(year == season, sub_field %in% field_names)
|
||||
|
||||
filtered_field_CI_data <- field_CI_data %>%
|
||||
filter(sub_field %in% field_names)
|
||||
|
||||
# CI <- map_df(field_names, ~ {
|
||||
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)
|
||||
|
||||
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) %>%
|
||||
mutate(DOY = seq(1, n(), 1),
|
||||
model = paste0("Data", season, " : ", field_names),
|
||||
season = season,
|
||||
sub_field = field_names)
|
||||
# }) #%>%
|
||||
#{if (length(field_names) > 0) message("Done!")}
|
||||
|
||||
return(CI)
|
||||
}
|
||||
|
||||
## Extracting the correct CI values
|
||||
# Data_2021 <- map(pivot_select_model_Data_2021, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2021)) %>% list_rbind() # nolint: line_length_linter.
|
||||
# message('2021')
|
||||
# Data_2022 <- map(pivot_select_model_Data_2022, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2022)) %>% list_rbind()
|
||||
# message('2022')
|
||||
# Data_2023 <- map(pivot_select_model_Data_2023, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2023)) %>% list_rbind()
|
||||
# message('2023')
|
||||
Data_2024 <- map(pivot_select_model_Data_2024, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2024)) %>% list_rbind()
|
||||
message('2024')
|
||||
|
||||
|
||||
#CI_all <- rbind(Data_2023, Data_2024)
|
||||
CI_all <- Data_2024
|
||||
message('CI_all created')
|
||||
#CI_all <- Data_2023
|
||||
|
||||
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')
|
||||
# nolint end
|
||||
print("All CI values extracted from latest image.")
|
||||
saveRDS(pivot_stats2, here(cumulative_CI_vals_dir,"combined_CI_data.rds")) #used to save the rest of the data into one file
|
||||
|
|
|
|||
141
r_app/interpolate_growth_model.R
Normal file
141
r_app/interpolate_growth_model.R
Normal file
|
|
@ -0,0 +1,141 @@
|
|||
|
||||
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 <- "sony"
|
||||
}
|
||||
|
||||
|
||||
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
|
||||
planet_tif_folder <- here(laravel_storage_dir, "merged_tif")
|
||||
merged_final <- here(laravel_storage_dir, "merged_final_tif")
|
||||
|
||||
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("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()
|
||||
|
||||
|
||||
# #2021
|
||||
# pivot_select_model_Data_2021 <- harvesting_data %>% filter(Year == 2021) %>% pull(field)
|
||||
#
|
||||
# pivot_select_model_Data_2022 <- harvesting_data %>% filter(Year == 2022) %>% pull(field)
|
||||
|
||||
# pivot_select_model_Data_2023 <- harvesting_data %>% filter(year == 2023) %>% filter(!is.na(season_start)) %>% pull(sub_field)
|
||||
|
||||
pivot_select_model_Data_2024 <- harvesting_data %>% filter(year == 2024)%>% filter(!is.na(season_start)) %>% pull(sub_field)
|
||||
|
||||
# pivots_dates_Data_2022 <- pivots_dates0 %>% filter(!is.na(season_end_2022))
|
||||
# pivot_select_model_Data_2022 <- unique(pivots_dates_Data_2022$pivot_quadrant )
|
||||
#
|
||||
# pivots_dates_Data_2023 <- pivots_dates0 %>% filter(!is.na(season_start_2023))
|
||||
# pivot_select_model_Data_2023 <- unique(pivots_dates_Data_2023$pivot_quadrant)
|
||||
#
|
||||
# pivots_dates_Data_2024 <- pivots_dates0 %>% filter(!is.na(season_start_2024))
|
||||
# pivot_select_model_Data_2024 <- unique(pivots_dates_Data_2024$pivot_quadrant)
|
||||
|
||||
extract_CI_data <- function(field_names, harvesting_data, field_CI_data, season) {
|
||||
# field_names = "4042902"
|
||||
# field_names = "1.13A"
|
||||
# harvesting_data = harvesting_data
|
||||
# field_CI_data = pivot_stats_long
|
||||
# season= 2024
|
||||
|
||||
filtered_harvesting_data <- harvesting_data %>%
|
||||
na.omit() %>%
|
||||
filter(year == season, sub_field %in% field_names)
|
||||
|
||||
filtered_field_CI_data <- field_CI_data %>%
|
||||
filter(sub_field %in% field_names)
|
||||
|
||||
# CI <- map_df(field_names, ~ {
|
||||
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)
|
||||
|
||||
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) %>%
|
||||
mutate(DOY = seq(1, n(), 1),
|
||||
model = paste0("Data", season, " : ", field_names),
|
||||
season = season,
|
||||
sub_field = field_names)
|
||||
# }) #%>%
|
||||
|
||||
return(CI)
|
||||
}
|
||||
|
||||
## Extracting the correct CI values
|
||||
# Data_2021 <- map(pivot_select_model_Data_2021, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2021)) %>% list_rbind()
|
||||
# message('2021')
|
||||
# Data_2022 <- map(pivot_select_model_Data_2022, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2022)) %>% list_rbind()
|
||||
# message('2022')
|
||||
# Data_2023 <- map(pivot_select_model_Data_2023, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2023)) %>% list_rbind()
|
||||
# message('2023')
|
||||
Data_2024 <- map(pivot_select_model_Data_2024, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2024)) %>% list_rbind()
|
||||
message('2024')
|
||||
|
||||
|
||||
#CI_all <- rbind(Data_2023, Data_2024)
|
||||
CI_all <- Data_2024
|
||||
message('CI_all created')
|
||||
#CI_all <- Data_2023
|
||||
|
||||
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')
|
||||
|
||||
|
|
@ -87,64 +87,10 @@ print(dates)
|
|||
#load pivot geojson
|
||||
# pivot_sf_q <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect()
|
||||
|
||||
raster_files <- list.files(planet_tif_folder,full.names = T, pattern = ".tif")
|
||||
|
||||
filtered_files <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>%
|
||||
vrt_files <- list.files(here(daily_vrt),full.names = T)
|
||||
vrt_list <- map(dates$days_filter, ~ vrt_files[grepl(pattern = .x, x = vrt_files)]) %>%
|
||||
compact() %>%
|
||||
flatten_chr()
|
||||
head(filtered_files)
|
||||
|
||||
# filtered_files <- raster_files #for first CI extraction
|
||||
|
||||
# create_mask_and_crop <- function(file, field_boundaries) {
|
||||
# message("starting ", file)
|
||||
# CI <- rast(file)
|
||||
# # names(CI) <- c("green","nir")
|
||||
# message("raster loaded")
|
||||
#
|
||||
# CI <- CI[[2]]/CI[[1]]-1
|
||||
# # CI <- CI$nir/CI$green-1
|
||||
# message("CI calculated")
|
||||
# CI <- terra::crop(CI, field_boundaries, mask = TRUE) #%>% CI_func()
|
||||
#
|
||||
# new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif"))
|
||||
# writeRaster(CI, new_file, overwrite = TRUE)
|
||||
#
|
||||
# vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
|
||||
# terra::vrt(new_file, vrt_file, overwrite = TRUE)
|
||||
#
|
||||
# return(CI)
|
||||
# }
|
||||
|
||||
|
||||
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")
|
||||
|
||||
vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt"))
|
||||
|
||||
if(emtpy_or_full[1,] > 10000){
|
||||
vrt_list[vrt_file] <- vrt_file
|
||||
|
||||
}else{
|
||||
file.remove(vrt_file)
|
||||
|
||||
}
|
||||
|
||||
message(file, " processed")
|
||||
gc()
|
||||
}
|
||||
|
||||
# list_global <- list_global %>% flatten_chr()
|
||||
vrt_list <- vrt_list %>% flatten_chr()
|
||||
|
||||
total_pix_area <- rast(vrt_list[1]) %>% terra::subset(1) %>% setValues(1) %>%
|
||||
crop(field_boundaries, mask = TRUE) %>%
|
||||
|
|
|
|||
Loading…
Reference in a new issue