From 1f194c567020a62cf334f2acadcdf9a48d60c266 Mon Sep 17 00:00:00 2001 From: Martin Folkerts Date: Mon, 1 Jul 2024 13:39:21 +0200 Subject: [PATCH] wip --- interpolate_growth_model.sh | 29 ++++++ r_app/ci_extraction.R | 155 ++++++++++++------------------- r_app/interpolate_growth_model.R | 141 ++++++++++++++++++++++++++++ r_app/mosaic_creation.R | 58 +----------- 4 files changed, 230 insertions(+), 153 deletions(-) create mode 100755 interpolate_growth_model.sh create mode 100644 r_app/interpolate_growth_model.R diff --git a/interpolate_growth_model.sh b/interpolate_growth_model.sh new file mode 100755 index 0000000..d694f4d --- /dev/null +++ b/interpolate_growth_model.sh @@ -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 \ No newline at end of file diff --git a/r_app/ci_extraction.R b/r_app/ci_extraction.R index 5695fa0..eac9719 100644 --- a/r_app/ci_extraction.R +++ b/r_app/ci_extraction.R @@ -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 diff --git a/r_app/interpolate_growth_model.R b/r_app/interpolate_growth_model.R new file mode 100644 index 0000000..dd6834e --- /dev/null +++ b/r_app/interpolate_growth_model.R @@ -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') + diff --git a/r_app/mosaic_creation.R b/r_app/mosaic_creation.R index 055604b..b132c97 100644 --- a/r_app/mosaic_creation.R +++ b/r_app/mosaic_creation.R @@ -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) %>%